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Mh, Abstract. The imbalanced superfluid state of spin-1/2 fermions with s-wave 

^ ' pairing is numerically studied by solving the Bogoliubov-de-Gcnncs equation at zero 

temperature in an annular disk geometry with narrow radial width. Two distinct 
types of systems are considered. The first case may be relevant to heavy fermion 
superconductors, where magnetic field causes spin imbalance via Zeeman interaction 
i-rt . and the system is studied in a grand canonical ensemble. As the magnetic field 

!Z| ' increases, the system is transformed from the uniform superfluid state to the Fuldc- 

Ferrell-Larkin-Ovchinnikov state, and finally to the spin polarized normal state. The 
second case may be relevant to cold fermionic systems, where the numbers of fermions 
of each species are fixed as in a canonical ensemble. In this case, the groundstatc 
^. . depends on the pairing strength. For weak pairing, the order parameter exhibits a 

^^ | periodic domain wall lattice pattern with a localized spin distribution at low spin 

imbalance, and a sinusoidally modulated pattern with extended spin distribution at 
high spin imbalance. For strong pairing, the phase separation between superfluid state 
and polarized normal state is found to be more preferable, while the increase of spin 
f ^ ■ imbalance simply changes the ratio between them. 
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1. Introduction 

In a conventional BCS theory, the normal state has a Fermi surface common to both 
spin-up and spin-down electrons and the Cooper pair has a zero total momentum. More 
than forty years ago, Fulde and Ferrell[T](FF), Larkin and Ovchinnikovp] (LO) proposed 
independently the pairing mechanism for the mismatched Fermi surfaces due to the spin 
imbalance. In the FF state, a spin up electron with momentum k is bounded with a spin 
down electron with momentum —k + q, thereby the Cooper pair has a net momentum 
q which is determined by the imbalance between two Fermi surfaces. Therefore the 
order parameter is characterized by a single momentum q, which can be written as 
A(f) = A Q e iq ' r with a uniform magnitude A . If considering the composition of two 
momenta, q and —q, one gets the LO state where the order parameter is real with its 
magnitude oscillating periodically in space. 

In condensed matter physics, the spin imbalance can be generated by applied 
magnetic fields. However the condition for the FFLO state to be observed is 
quite stringent on the superconducting materials. Roughly speaking, there are three 
requirements (i) low T c , so that the magnetic field needed to imbalance the spin 
population is accessible; (ii) the orbital effect of magnetic field is weak enough to 
avoid pair breaking before the Zeeman splitting takes effect; (iii) clean limit, i.e., the 
mean free path of electron should be much longer than the correlation length, since the 
FFLO state is easily destroyed by impurities. Some of heavy fermion superconductors 
are good candidates to fulfill these requirements (for a review see Ref. [3]). There 
was recent indications that CeCoIn 5 indeed exhibits the FFLO state [U [3 [8] [9]. 
That compound is a quasi-two-dimensional heavy fermion superconductor with a d- 
wave pairing. In the cold fermionic atom system with different hyperfine spins, the 
spin population imbalance between different hyperfine spins can be easily controlled 
by applying radio frequency field. Recently the imbalanced superfluid state has been 
realized in these cold neutral atom systems [TUJ [HI H21 [131 EH] > an d the possible spatially 
modulated superfluid phases in these systems are studied in Ref. [15], [16]. It is noted that 
the particle number of different species may be controlled directly in systems like cold 
atoms, and in superconductors the spin imbalance is generated by the external magnetic 
fields, which may correspond to two different thermodynamic conditions, respectively. 

In a recent theoretical study [IT], it was found that in the harmonically trapped 
polarized fermionic atoms in a two-dimensional (2D) optical lattice, the insulating core 
is surrounded by a superfluid shell at high atom densities with pairing parameter 
modulated in the circumferential direction. Since some of important physics may be 
explained by the quasi-one-dimensional (quasi-lD) shell, it is thus interesting to study 
further the FFLO with more details in a quasi- ID system. The possible angular FFLO 
state in a toroidal trap has also been investigated in a very recent study [18]. In the 
present paper, we consider a quasi-lD annular disk with narrow enough radial width, 
so that the radial modulation of the order parameter might result in a quite large radial 
gradient of order parameter which increases the system energy considerably according 
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to the Ginzburg Landau(GL) theory. Therefore the oscillation of pairing amplitude is 
suppressed in radial direction, and restricted only in circumferential direction. In a 
large 2D system, the order parameter oscillation has more freedom and can happen in 
arbitrary directions. In the presence of inhomogeneity the modulation direction may 
vary in space which leads to irregular pattern of order parameter. Therefore it may be 
easier to observe regular oscillations of the pairing amplitude in a quasi- ID system than 
in the 2D film. 

In this paper, we consider two distinct systems. The first one may be relevant 
to heavy fermion superconductors, where the electrons spins interact with an external 
magnetic field via the Zeeman coupling. The second system may be related to the 
cold fermionic atoms, where the number of fermions of each spin is fixed. We employ a 
grand canonical ensemble to study the first system and a canonical ensemble to study the 
second system. We solve the Bogoliubov-de-Gennes (BdG) equation at zero temperature 
numerically for the above quasi-lD systems. Our main results can be summarized 
below. In the first case, as the magnetic field increases, the ground state is transformed 
from a uniform superfluid state to the sinusoidally modulated LO state, and then to 
a spin polarized normal state. In the second case, the ground state depends on the 
pairing strength. For weak interactions, the order parameter exhibits a periodic domain 
wall lattice pattern with a localized spin distribution for low spin imbalance, and a 
sinusoidally modulated pattern with extended spin distribution for high spin imbalance. 
For strong interactions, the phase separation between superfluid state and polarized 
normal state is found to be more preferable, while increase of spin imbalance simply 
extends the spatial region of the normal state. The paper is organized as follows. In 
Sec. II, we study the exact ID case. In Sec. Ill, we present our results for annular disk 
geometry. The conclusion is given in Sec. IV. 

2. Imbalanced Superfluid State in One-dimensional Ring 

2.1. One-dimensional BdG Equation 

Before exploring the properties of imbalanced superfluid in the annular disk geometry, 
we first consider the ID ring, which may be viewed as the limiting case where the 
disk width is so narrow that only one radial mode is relevant. This case has been 
studied by a number of authors. In the mean field(MF) level, a rigorous analysis for 
the ID BdG equation is given in Ref. [19J in the presence of a magnetic field. In terms 
of ID Luttinger liquid theory the imbalanced superconducting state is also elucidated 
by Yang [20], and very recently, the density matrix renormalization group algorithms 
are implemented on the ID negative-?/ Hubbard model to explore the FFLO state in 
Refs. [2H |22l [23J [21] . The cold fermionic gases with attractive interaction and population 
imbalance are studied theoretically in Ref. [25] and and Ref. |26j. 

In this subsection, we follow the MF treatment to give a brief description to the 
ID imbalanced superfluid state. We consider a canonical ensemble and fix the number 
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of fermions of different species. Although only the quasi-long range order may exist 
in ID system, the MF approach presented in this section is helpful to understand the 
imbalanced superfluid in 2D annular disk geometry shown in later sections. 
The mean field Hamiltonian for a ID interacting system reads 



H = J dx(£#U*) (^^) U*) 



A(x)'ipUx)ijjUx) + h.c. 



~^2^a[ dx^ a ( X )^a( X ) - N a 



A(x) = g (faWhix)) . (1) 

i/j a { x ) is the fermion annihilation field at position x with spin index a, A(x) is the 
fermion pairing field, m is the mass of the particle, and g < is the attractive interaction 
strength. fi a are the Lagrangian multipliers or the chemical potentials, which are used 
to fix the numbers of fermions of different spins at N^ and iVj_, respectively. 

Eq. (CQ) has the similar form to the well known Su-Schrieffer-Heeger(SSH) model 
[27] for polyacetylene, which describes a ID electron system coupled to phonons. In 
this system, when the phonon fields are condensed in opposite phases at the two ends 
of the ID string, there are possible soliton excitations with zero energy in the fermion 
spectrum. The soliton excitations are also possible in the ID superfluid Hamiltonian 
Eq. (CQ), where the MF pairing parameter A(x) can mimic the phonon field in the SSH 
model, which is shown briefly below. More details can be found, e.g., in Ref. [19]. For 
simplicity we take /^ — jUj, = fi, which determines the Fermi momentum kp = y/2mfj,/h. 
The low energy physics is described by quasiparticles around the two Fermi points ±fcp, 
i.e., the following decomposition is allowed 

Mx) ~ e ik » x R ff {x) + e- tkFX L a (x) (2) 

with left and right movers defined as 



Ra{x) = ^ i>v{k + k F ) r - 

-A<fc<A V L 

La(x) = ^ 4>a{k-k F )— . (3) 



-A<fc<A v 

A is a suitable momentum cutoff. These quasiparticle operators satisfy the standard 
anti-commutation relations, i.e., 



j (j,(j i 



{R a ,Rl} = {L CT ,Ll l } = 6 

and all the other anti-commutators are zero. Substituting Eq. [2] into Eq. [TJ and 
neglecting the fast oscillation terms (oc exp(±2ikpx)), one obtains the following two 
Hamiltonians to the linear order of fe, 



Hi = Kvf / dx : Ri{—id x )R-\ : — : L,U—id x )Ln 
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+ dxA{x){R\L\+L l Ri i ) 
H 2 = hv F / dx : R\(—id x )R± : — : Ll(—id x )L^ : 

+ J dxA{x){l\R\ + R X L^ . (4) 

Here : A : denotes normal ordering of A and Vp means the positive Fermi velocity. In 
the following Hvf is taken as unit. Hi and H2 are commutative with each other, and 
connected through the gap equation 

A(x) = g (R^ + L t Ri) . (5) 

The order parameter A(x) is assumed to be real. Eq. [5] shows that the pairing takes 
place either between Lf and i2j_, or between Lj_ and Rp Actually, (R^L^) = (L^R-^) by 
symmetry. Formally, one may have H ~ Hi + H 2 — f dx\A(x)\ 2 /g, but it is emphasized 
that Hi 2 only describe the low energy excitations near the Fermi surface. 

Let's consider only Hi with a twisted A(x), i.e., A(— 00) = — A(oo) = A . As 
shown by Jackiw and Rebbi[28], there is at least one zero mode 70^ in the middle of the 
gap, which is localized in space and reads 



Tot oc 



dxF(x)[R^(x) -ilUx)] 



F(x) oc exp 



.1; 



dx'A(x') 



(6) 



It is easy to verify the commutation relation [Tot , -Hi] — 0. Besides this localized zero 
mode, we also have other quasiparticle excitations Tna i n the continuum region, where n 
and a are the energy level and spin indices, respectively. Assuming all of them constitute 
a complete representation of the Hamiltonian Hi, the lowest energy states are doubly 
degenerate in the presence of an order parameter with kink pattern, which is the spinless 
vacuum of the quasiparticles Tna together with the zero mode Tot being either filled or 
empty. Similar analysis is also valid for the H 2 branch, for which one can find that the 
zero mode has the form 

Tot oc / dxF(x)[Ri(x) +iLl(x)] 

which satisfies [701, #2] = 0. 

In terms of R a and L a , the total particle number A^ and total spin operator S can 
be written as 

N = A> T + N h S = N^-N l 

N a = f dx[. RlRc : + : LlL. :] , 

where the fast oscillating terms are neglected. Note that the quasiparticle operators R a 
and L a can only describe the low energy physics, hence the operator A CT with normal 
ordering only measures the particle number relative to the Fermi surface. Obviously, 
unlike the SSH model[27j and the Jackiw- Rebbi model[28], the charge conservation is 



Imbalanced superfluid state in an annular disk 



6 



broken in the BCS theory, therefore one can not tell how many charges the soliton can 
carry Despite this fact, the total spin is still a conserved quantity in our MF treatment, 
therefore each zero mode may carry half spin as an analog to the half charge investigated 
in Ref . [271 [28] . But in practice only one spin can be observed at the kink of A(x), since 
there are two branches of fermions(iJi and H2). To observe the half spin, one must 
get rid of the fermion doubling problem. Nevertheless, this provides a mechanism to 
accommodate excess spins with zero energy. The total energy of the soliton measured 
relative to the uniform BCS state is computed to be 2A /7r[29j [30], which is less than 
the superfluid gap. 



2.2. From Soliton Lattice-like LO state to sinus oidally-varying LO State 

For equally populated species N^ = iVj_, the lowest energy state is obviously the BCS 
state with uniform pairing gap. If one spin is flipped from downward to upward, i.e., 
iVf + 1 up spin and iVj_ — 1 down spin, a pair of soliton and anti-soliton is developed 
to store these two excess spins. We define the spin imbalance n to be (JVf — iVjJ/2 for 
spin 1/2 particle. A typical soliton and anti-soliton pair is plotted in Fig. (Ilk), which 
is obtained by numerically solving Eq. (TjQ) in a ring, where we use the angle 9 = 2irx/L 
as the coordinate. Due to the periodic boundary condition, a single soliton can not 
exist freely so that it must co-exist with an anti-soliton as a pair with the same width 
£. We call these soliton states with each spin per soliton (anti-soliton) as ideal soliton 
state. Note that since the order parameter is real, this state is also a kind of LO state. 
Actually, all the self-consistent solutions shown in this paper have real order parameters 
which minimize the energy, and therefore they are LO state. In the following sections, 
we omit "LO" for the sake of brevity. 
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Figure 1. Angle distribution of pairing order parameter in an ideal soliton state. The 
order parameter is measured in unit of Ao which is the value of order parameter in the 
uniform state. From top to bottom, total spin imbalance is 1, 6, and 14. Open symbols: 
numerical results; solid lines: fitting function Atanh(cosn#/£) with two parameters A 
and £. 
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With the increase of the flipped spins, more soliton and anti-soliton pairs are 
generated. Thus we get the soliton lattice state with n pairs of soliton and anti-soliton 
as long as the system is in the dilute limit by which we mean n£ <C 2n, here the 
soliton width £ is measured in unit of the angle. In the dilute limit, the solitons 
are well separated from each other, which has two consequences (i) all the midgap 
states have zero energy, and (ii) each soliton or anti-soliton carries exactly one localized 
spin. According to these two properties, we distinguish soliton lattice state from the 
sinusoidally modulated state, where the spin imbalance n is too large to satisfy n£ < 2tt 
and solitons overlap considerably with each other. Then the energy spectrum of the 
midgap states has a dispersion described by the Bloch theorem for a periodic lattice. 
Such a scenario from soliton lattice to sinusoidally varying state has also been addressed 
in Ref. [31 j from the viewpoint of GL theory. The pairing parameter for both states 
can be described perfectly by the fitting function Atanh(cosn^/£) L0 with A and £ to 
be determined, which is shown in Fig. HJ 

We now introduce two spin distribution functions, local spin distribution Sl(8) = 
\(ilr^{9)%[)}{9) — ■ijj^(9) , ijji(6)) as well as integrated spin distribution Si(9), 



Si(0) 



S L (6')d6' . 



(7) 



As shown in Fig. [2J the localization of spin density in the soliton lattice state manifests 
itself in the plateau features of the function Si{6). For the sinusoidally modulated state, 
the plateaus disappear due to the derealization of spins. 
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Figure 2. Spin distribution in soliton lattice state for the system with spin imbalance 
6 (a) and 14 (b). Solid lines: local spin distribution; Dash lines: integrated spin 
distribution. Inset of (b) shows a zoomed figure around a plateau. 



% The soliton lattice pattern of pairing parameter can be described by Jacobi elliptic function as done 
in Ref. |19j , but we do not take that expression for the sake of simplicity. 
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2.3. Deformed Soliton 

Here we introduce Q to denote the number of spins per soliton/antisoliton. In the 
previous subsections, we focused on the state with only one spin(Q = 1) per soliton. 
Now we study the case for Q > 2, which we call deformed soliton state. Firstly, let us 
consider the case for odd Q. The order parameter of a deformed soliton state for Q = 3 
is plotted in Fig. ^ solid lines), which corresponds to 6 excess spins in total. Note that 
these 6 spins can also be stored in 3 ideal soliton-antisoliton pairs(dashed lines). Hence, 
we need to compare their energies numerically. It turns out that the deformed soliton is 
energetically favorable for strong interaction, while the ideal soliton state is preferable 
for weak interaction. Note that the deformed soliton found in this article has Q nodes 
in a narrow region. In fact Q spins can also be accommodated by a special soliton with 
only one nodes, which is described by A tanh(x/£) with A £ = (<2 + l)/2(see Ref.[30]). 
however one can show that this solution is not energetically favored by comparing its 
energy and that of the corresponding well separated multi-soliton state. 

In Fig. EH the upper panel corresponds to a strong interaction case where the three 
spins are squeezed in a very narrow region with width comparable to that of an ideal 
soliton £. The total width is then estimated to be around 2£, which is much smaller 
than the width 6£ for the ideal soliton state. Thus, one can reasonably believe that the 
deformed soliton state has lower energy. If the interaction strength g becomes weaker, 
as shown in the lower panel of Fig. [3], the deformed soliton with Q = 3 will inflate 
and its pattern is getting close to three ideal solitons. When g becomes weak enough, 
the deformed soliton can not be stable, and is transmuted into an ideal soliton lattice 
state. The pairing order parameter shown in Fig. [3] can be perfectly fitted with function 
A [tanh(cos(6> — #o)/£o) — tanh(cos(6>)/£ + tanh(cos(6 l + #o)/£o] with three parameters 
A , | and 9 . 




0.4 0.6 

0/(2tt) 



Figure 3. Pairing parameter in deformed soliton with Q = 3 (solid lines) and ideal 
soliton with Q = 1 (dashed lines). Upper (lower) panel corresponds to the strong 
(weak) pairing interaction g. 



Note that the order parameter has a sign change (— l) s after crossing Q ideal 
solitons and antisolitons. Therefore, if Q is odd, a deformed soliton can be continuously 
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transmuted into Q ideal solitons, but this is not true for even Q due to the mismatched 
boundary condition of A(x). In addition, the energy of a deformed soliton with even Q 
is not energetically favorable in our numerical calculations. Therefore, we do not need 
to consider the case for even Q. 

2.4. Effect of Magnetic Field 

So far we only consider the system with fixed particle number, and have not included 
the magnetic field in our analysis. Since the total spin is a good quantum number, 
the effect of magnetic field can be easily estimated by simply adding Zeeman energy 
—fiBh(N^ — iV|). Obviously, the state with more excess spins gains magnetic energy, 
however it is at the cost of the deformation of pairing gap which loses the condensation 
energy. Therefore, the ground state should correspond to an optimized value of spin 
imbalance. 

Let n = (iVj- — iVjJ/2 be the spin imbalance, and the corresponding ground state 
energy be denoted by E(n). The energy of the BCS state without spin imbalance is thus 
E(0). Given an external magnetic field h, we then need to find the lowest free energy 
for all possible n's, i.e., minimize E(n) — 2n[iBh with respect to n, which leads to an 
optimal spin imbalance n c . 

To this purpose, we define the energy cost per spin as 

e(n) = [E(n)-E(0)]/(2n), (8) 

which can also be regarded as the energy cost for creating one soliton. The numerical 
data of e(n) is plotted in Fig. HI As n increases, the adjacent kinks become closer, 
which enhances the hopping amplitude of spins between kinks and consequently favors 
the kinetic energy of spin transfer. However, at the same time, the pairing gap gets 
smaller, which reduces the condensation energy. Thus the interplay between these two 
mechanisms leads to the nontrivial pattern of e(n) in Fig. HI 




5 10 15 20 25 30 35 40 45 50 
spin imbalance 

Figure 4. Average energy per spin e(n) in Eq. ([8]) as a function of spin imbalance n. 
The dashed line is the first critical magnetic field h\. 



There is a critical value hi of the magnetic field, below which the magnetic energy 
can not support an ideal soliton, and the system remains in the uniform state. When 
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h > hi, the sinusoidally varying state with modulation frequency n c will become 
energetically favorable. n c can be determined by the minimum of 2ne(n) — 2np,Bh, 
alternatively, the optimal spin imbalance n c should satisfy 

<9(2n(e(n) — /Ab/i)) 



dn 



= (9) 

n=n c 

After a little algebraic analysis of Eq. (|9"1), one can see that n c increases as h increases. 
The first appeared n c is determined by e(n c ) = pshi which is far from zero as shown in 
Fig. H] and corresponds to a sinusoidally modulated state. 

3. Imbalanced Superfluid State in Annular Disk 

In this section we present our numerical results for imbalanced superfluid state in narrow 
annular disk with inner radius Ri and outer radius R 2 . The radial width R 2 — Ri is small 
enough to avoid the modulation of order parameter along the radial direction. In the 
numerical calculation, we use the ratio p = (R 2 — Ri) / Ri to characterize the geometry of 
annular disk. Since g has the dimension of [energy] -[length] 2 , a dimensionless quantity 
g = g/(ir(R% — -R 2 )/i) is introduced to represent the interaction strength. The BdG 
equation is solved in momentum space. Most of the results in this section are based 
upon the diagonalization of Hamiltonian in a Hilbert space with dimensionality 3500 
and 11 radial modes involved. 

3. 1 . Fixing Particle Number JVj and N^ 

3.1.1. Ideal Domain Wall For small spin imbalance, one should get domain walls 
as an analog of solitons in ID case, and the excess spins are attached to the domain 
walls. It is natural to ask what is the optimal number(Q) of spins per domain wall. 
To answer this question, we first consider an ideal geometry, i.e., a narrow strip with 
periodic boundary condition in both x and y directions, but with length L x ^> L y . 
This simplified model reads 

/ "-2 



H = dxdy[ipl I fi a J i> a 

9 

A(x,y)=g^ l ^ . (10) 

The ideal domain wall pattern of A(x,y) is independent of y, and has the form 
A(x,y) = A tanh(a;/£o) which implies the pairing momenta in y direction are always 
q and —q. The Hamiltonian in Eq. ( TTUl) can be divided into many ID branches with 
respect to the discrete momenta q = (2ir/L y ) x integer in y direction, 



H q ~ dx^l^x) (ZjL-Htf) 4,T<X> 
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+ il>l 



■q,l 



X 



2^ " ^ 



$-q,l(x) 



+ A(i)C,(j;)+AV-,,i(x)it(*)] • 



:n) 



Note that A(x) is contributed from all ID branches, and the g-dependent chemical 
potential reads \i qa — // a — (hq) 2 /(2m), which are determined by the particle numbers 
N a . Each g-mode with \x q > can accommodate one spin per soliton. Therefore, we 
can estimate the optimal spin filling Q of each ideal domain wall to be the number 
of g-modes buried under the FS. The optimal filling for the annular disk with open 
boundary condition in the radial direction can also be estimated similarly by counting 
the number of energy modes under the FS. 

Similar to the ID ring, one expects a crossover from an ideal domain wall like LO 
state to the sinusoidally-varying LO state with increasing spin imbalance in the weak 
interaction case. Since A(r, 9) now depends on r, we plot the angle dependence of A(r, 9) 
at radius r = (R\ + R-i)j2 in Fig. [5l The full spatial dependence of A(r, 9) is plotted in 
2D contour in Fig. [7J where one can find its radial dependence is nearly uniform. The 
spin density s(r,9) is also a function of r and 9. By integrating s(r,9) over r, we can 
define angle dependent local spin distribution Sl(9), and angle dependent integrated 
spin distribution Sj(9), as following, 

S L (9) = / rdrs(r, 9) 



R\ 



Sj(9) 



S L (9)d9' . 



(12) 



Sl and Si are plotted as functions of 9 in Fig. EJ which shows clearly that the spin 
distribution are localized in the domain wall state, and delocalized in the sinusoidally- 
varying LO state. 
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Figure 5. Angle dependence of pairing order parameter at radius (i?i + R,2)/2. (a): 
Domain wall lattice state with total spin 28; (b): sinusoidally-varying LO state with 
total spin 70. The optimal filling per domain wall is Q = 7. The system parameter 
,9 = 0.4. 
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Figure 6. Spin distribution in domain wall lattice state with spin imbalance 28 (upper 
panel) and in sinusoidally- varying state with spin imbalance 70 (lower panel). Dashed 
lines: integrated spin distribution Si (9), and solid lines: local spin distribution Sl{9)- 
The system parameters are the same as in Fig. [5) 
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Figure 7. (Color online.) Contour plot of order parameter. The excess spin equals to 
28 and the optimal filling in this case is Q = 7, hence four pairs of domain walls are 
needed to store these excess spins, the system parameters are the same as in Fig. [5j 



3.1.2. Deformed Domain Wall and Phase Separation As in the ID ring, we also 
encounter the deformed domain wall state, for which there can be more spins than 
the optimal filling Q squeezed in one domain wall. These deformed domain wall states 
are stabilized by the strong pairing interaction. We plot the order parameter A(9, r) and 
local spin distribution Sl(0) in Fig. [8], which shows that when the spin number exceeds 
the optimal filling, instead of creating more ideal domain walls, the spin polarized 
regions are simply enlarged. Note that in the polarized region there is still a small 
pairing oscillation like a mini sinusoidally- varying LO state in order to further lower the 
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potential energy. These deformed domain wall states (see Fig. |Sb) are then considered as 
a kind of phase separation state, where the polarized normal state with small fluctuating 
order parameter is separated with the fully pairing phase without spin imbalance. 
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Figure 8. Deformed domain wall [(a) and (b)] and phase separation [(c) and (d)] 
solutions. We plot the order parameter in (a) and (c), and spin distribution in (b) and 
(d). The spin imbalance is 21 for (a) and (b), and 77 for (c) and (d). The optimal spin 
filling Q — 7. The interaction strength is g ~ 6.9 x 10~ 4 . 



3.1.3. Quasiparticle Density of States We compute the quasiparticle density of states 
(DOS) in this section which can describe the low energy excitations of various ground 
states. In our calculation the Zeeman energy is not included, which corresponds to 
the situation with fixed particle numbers. We find that, for the domain wall lattice 
state there is a zero energy peak in the quasiparticle DOS. As the spin imbalance is 
increasing, the number of domain walls grows and it results in the enhancement of 
the zero energy peak. These zero modes can also be understood from the aspect of 
Andreev reflection [32], since the 7r-phase difference between two superfluids allows an 
Andreev bound state located at the domain walls. In the phase separation case, the 
system mimics a superconductor-normal metal-superconductor junction. By increasing 
the width of normal metal region, more Andreev resonance states enter into the gap 
with nonzero energy. These energy levels then distribute evenly in the gap, which form 
a flat quasiparticle DOS in the superconducting gap. 

The above theoretical analysis is in good agreement with the numerical results 
presented in Fig. [9] The DOS of BCS state is zero in the gap. When increasing the 
spin imbalance in the ideal domain wall lattice state, the peak of DOS centered around 
zero becomes higher, which means more domain walls are created. Whereas in the case 
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of phase separation, the DOS in the gap is quite flat due to the presence of polarized 
normal state. 
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Figure 9. (Color online.) Quasiparticle density of state for different ground states. 
The Zeeman energy is not included in this figure. The red solid line is for the uniform 
BCS state, the green long dashed line and the blue short dashed line are for the domain 
wall lattice states, the dotted pink line is for the sinusoidally modulated LO state, and 
the cyan dot-dashed line corresponds to the phase separation state. 



3.2. Fixing Chemical Potentials /x-f and /xj_ 

In this subsection, we show the numerical results in the grand canonical ensemble with 
fixed chemical potentials. For weak magnetic field(2/i#/i = //-j- — fi[), the Zeeman energy 
is not enough to break the s-wave Cooper pairs, so the system retains the uniform BCS 
state. Until the magnetic field h exceeds its first critical value hi, the sinusoidally- 
varying LO state emerges. As the magnetic field is further increased, the modulation 
frequency of the order parameter becomes larger while its magnitude is reduced, until 
the system enters into the normal state at the second critical magnetic field hi- We plot 
modulation frequency as a function of h in Fig. [10l where one can find plateaus, since 
there should be integral pairs of domain walls in a ring geometry. 

The phase separation(deformed domain wall) state can not be a ground state 
in the homogeneous magnetic field, except at the critical value hi of magnetic field. 
Furthermore, unlike the case of fixing particle number, there is no continuous crossover 
from domain wall state to the sinusoidally-varying state. The onset frequency at the 
critical magnetic field hi is finite and large enough to form a sinusoidally-varying LO 
state. The reason is that, to sustain a single domain wall, its magnetic energy gain 
must fully compensate the energy loss due to the deformation of pairing gap. In 
such a case there can be more domain walls. However the overlap of domain walls 
suppresses the pairing gap inevitably, which causes the loss of the condensate energy(see 
sec. 12.41) . At the balance point of these two processes, sinusoidally-varying state shows 
up accompanied with delocalized spins. 
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Figure 10. Frequency of pairing modulation as a function of magnetic field. We set 
p = 0.2, and g = 5 x 10~ 4 . [Ib is the Bohr magneton, and g-factor of electron is taken 
as 2. 



4. Conclusion 



We have investigated the imbalanced superfluid state in annular disks and ID rings by 
solving the BdG equation in the momentum space at zero temperature. A key issue of 
imbalance superfluid is how to accommodate the excess spins by adjusting the pairing 
gap A(r). There are several possibilities, e.g. the LO state with periodically oscillated 
order parameter and the phase separation state. We show that these states are stable 
under different conditions. 

Firstly, we have studied the case with fixed fermion numbers, which may be relevant 
to cold atom systems. For low spin imbalance (still larger than the optimal spin filling Q 
per domain wall), the solitons in ID and domain walls in 2D are the ground states. The 
number of spins localized at each soliton or domain wall is quantized. When increasing 
spin imbalance, more and more domain walls(solitons) occur and overlap with each 
other, and the sinusoidally-varying state emerges with delocalized spins. These two 
states are distinguished in this paper due to their different spin distribution. There 
should be a crossover between them if one tunes the spin imbalance continuously. The 
above argument is valid for weak interactions, whereas for strong interactions, the phase 
separation is the possible ground state, in which only the area of normal polarized state 
varies with the spin imbalance. This may serve as a criteria to distinguish the phase 
separation state and the periodically oscillating LO state. 

Secondly, we have addressed the case of fixing chemical potential fi and magnetic 
field h, which may be relevant to heavy fermion superconductors interacting with an 
external magnetic field via the Zeeman term. There are two critical magnetic fields hi 
and hi, which correspond to the transition from uniform BCS state to the sinusoidally- 
varying state, and from the sinusoidally-varying state to the normal state, respectively. 
It is stressed that the modulation frequency of pairing gap at h\ is quite large and the 
spin is delocalized, which characterizes a typical sinusoidally-varying state. 
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